Cite this Article
Wu Qi-Lin, Hou Mei-Ying, Yang Lei, Wang Wei, Yang Guang-Hui, Tao Ke-Wei, Chen Liang-Wen, Zhang Sheng. Parametric study of the clustering transition in vibration driven granular gas system. Chinese Physics B, 2020, 29(5): 054502  
Permissions
Parametric study of the clustering transition in vibration driven granular gas system
Wu Qi-Lin1, 2, Hou Mei-Ying3, Yang Lei1, Wang Wei1, 2, Yang Guang-Hui1, Tao Ke-Wei1, Chen Liang-Wen1, Zhang Sheng1, †
Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou 730000, China
University of Chinese Academy of Sciences, Beijing 100049, China
Key Laboratory of Soft Matter Physics, Beijing National Laboratory for Condense Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China

 

† Corresponding author. E-mail: zhangsheng@impcas.ac.cn

Abstract

A parametric study of the clustering transition of a vibration-driven granular gas system is performed by simulation. The parameters studied include the global volume fraction of the system, the size of the system, the friction coefficient, and the restitution coefficient among particles and among particle–walls. The periodic boundary and fixed boundary of sidewalls are also checked in the simulation. The simulation results provide us the necessary “heating” time for the system to reach steady state, and the friction term needed to be included in the “cooling” time. A gas-cluster phase diagram obtained through Kolmogorov–Smirnov (K–S) test analysis using similar experimental parameters is given. The influence of the parameters to the transition is then investigated in simulations. This simulation investigation helps us gain understanding which otherwise cannot be obtained by experiment alone, and makes suggestions on the determination of parameters to be chosen in experiments.

PACS: ;45.70.-n;;45.70.Mg;;81.30.Dz;
Keyword:;granular gas;;clustering;;phase diagram;
1. Introduction

We regard a system formed by a large number of macroscopic particles with a size larger than micrometers as a granular system. Granular matter is ubiquitous in our daily life. It is the most existing material form on earth except water. Examples such as sand, industrial products, food, grains, etc., are all considered as granular materials. Some natural phenomena such as sandstorm, avalanche, earthquake, debris flow, and other geological disasters are closely related to granular flow. The cognition of these phenomena and their mechanisms has been one of the foci of physicists’ interest.[1]

In the granular system, thermal fluctuation can be neglected when compared to particle kinetic energy. Thus, the granular system is considered as a system far from thermal equilibrium. It is an amorphous, discrete system with intrinsic energy dissipation. According to the volume fraction of the system and the kinetic energy of the particles, we can generally divide granular systems into “granular gas”, “granular liquid”, and “granular solid”. Unlike molecular gases, the particles in granular gas collide with each other in an inelastic way, and sedimentation occurs under gravity. Inelastic collision causes cluster formation in granular gas. The mechanisms of clustering have been of interest to physicists in recent years.[1–3] Understanding the mechanisms of the instability and the spontaneous formation of clusters in granular gas may help physicists establishing non-equilibrium statistical model.

To study the granular gas properties experimentally, energy input is necessary to balance the intrinsic energy dissipation due to the particle–particle inelastic collisions. The most commonly used method is to agitate the system mechanically by oscillation, that is, particles gain momentum through colliding with the oscillating-wall, the so-called “boundary heating” method. The agitated particles will then reach a “thermal” equilibrium through particle–particle collisions. However, sedimentation under gravity makes experimental investigation difficult. Therefore, it is necessary to conduct experimental study of the granular gas behavior either under microgravity[4,5] or by computer simulation.[6,7]

In a pioneering experiment, cluster forming in granular gas was reported by Falcon et al.[8] The system was mechanically driven into a steady-state by shaker to balance the internal energy dissipation. Clustering occurred when the global volume fraction ϕg exceeded a critical value, resulting in a decrease in both pressure and kinetic temperature.[9] The location of the cluster was far from the driving piston and ‘preferred’ to aggregate near the sidewalls.[10–12]

In event-driven simulation by Miller and Luding,[7] the growth of cluster was characterized by the power law decay of kinetic energy of the system. To identify a cluster, a statistical method Kolmogorov–Smirnov (K–S) test was proposed[13] to evaluate the deviation of the spatial distribution of particle from a uniform distribution. In the numerical work of Opsomer et al.,[14] the local Voronoi volume fraction ϕlocal exceeding 0.285 was considered as a criterion for cluster to appear. There seems to be consistency between the K–S test and criterion by ϕlocal.

In this work, we focus on the simulation investigation of the gas-cluster transition in piston driven granular gas system. The “heat” transport from oscillating-walls to the granular gas, and the “cooling” of the kinetic energy due to inelastic collisions among particles to cluster state are investigated. The K–S test and ϕlocal criteria used to establish a phase diagram are checked. We then investigate the influencing parameters of the transition curve in the phase diagram.

2. Methodology
2.1. Simulation model

The simulation is based on discrete element method (DEM).[15] We use Hertz nonlinear contact model to calculate the force between particles. The interaction between them is considered if and only if two particles collide. Taking the rotation into consideration, the particles’ equations of motion are as follows:

where mi, Ii, ri, and ωi are the mass, the moment of inertia, the displacement vector, and the angular velocity of particle i, respectively. is the contact force that the particle i is applied by the particle j. and are other forces and torques applied to the particle i, such as the gravity force. Generally, a normal force and a tangential force can compose

In the normal direction, the forces between particles include Hertz elastic force and damping force

where λn = (Ri+ Rj) – rij, rij = ri–rj is the overlap between i and j, Ri is the radius of particle i. The spring stiffness and the dissipative coefficient , where ε, Y*, m*, and R are the coefficient of restitution, equivalent elastic modulus, equivalent mass, and equivalent radius, respectively. And vnij is the normal component of the relative velocity vnij = vi – vj. In this work, the coefficient of restitution is a constant, independent of the relative velocity of two particles during the collision.

The tangential contact force is determined by the tangential component of the relative velocity and the normal contact force

where is tangential relative velocity, kt is tangent elastic coefficient, γt is tangent damping coefficient, utij is the tangential displacement between particles, which is related to the contact time τc of particles, utij = 0 at the first contact, and then its value can be obtained by

In addition, the yield condition is ensured by truncating the tangent displacement utij, and μ is the coefficient of friction.

2.2. Simulation system

The simulation system contains monodisperse spheres with a diameter of 1 mm. Other related parameters are listed in Table 1. The system contains a rectangular sample cell with two pistons moving in the x-direction to agitate the spheres in the cell, as shown in Fig. 1(a). We keep the Young’s modulus and Poisson ratio of the cell walls and particles identical, but keep the coefficient of sliding friction and the restitution coefficient variable (see Table 2). In the simulation, the pistons are in harmonic oscillations with constant amplitude and frequency. The two pistons are kept in opposite phases to avoid resulting in a coherent ‘bouncing’ state.[14] The total cell volume is V = Lx × 30 × 30 mm3, where Lx is the distance between the two moving pistons. The total number of particles in the cell is denoted as N. The simulation is done either with fixed boundary or with periodic boundary (PBC) in the y–z direction (see Table 2).

Fig. 1. (a) A schematic diagram of the simulation system. The volume of the system is V = Lx × 30 × 30 mm3, in which two pistons oscillate in amplitude A and frequency f. Different colors are used to indicate the Voronoi volume fraction ϕlocal of each particle. (b) A slice (1 mm thickness) chosen in (a).
Table 1.

The parameters used in simulation.

.
Table 2.

Seven sets of simulation parameters, where μpp, μpw, ε are respectively the coefficient of friction between spheres, the coefficient of friction between spheres and walls, and the coefficient of restitution. No. 1: reference; No. 7: freely cooling.

.

The maximum mean distance between two particles in the cell is given by . The initial condition of the particles is prepared nearly uniform in spatial distribution with the volume fraction ϕg = 4Nπr3/3V, where r is the radius of the particle. The initial velocity of each particle is set to be equal in magnitude, 0.01 m/s at t = 0, but the direction is random.

2.3. Criteria

One method to set the criterion of cluster from gas is by the two-sample K–S test.[16] In the method, the maximum difference between a cumulative distribution function F(x) of counts along the x-axis and the cumulative distribution function of a uniform function U(x) is defined as

The D is set to compare with a statistical threshold . The distribution is considered uniform if , otherwise it is non-uniform. In our case, we calculate cumulative distribution function F(x) of the number density distribution ρ(x) along the x-axis, and n is the number of sampling, n = Lx/dc (dc is the thickness of one slice, see Fig. 1(b)). Kα is a contrast threshold number.[17] The confidence level (1–α) chosen in this work is 0.95. When , we consider that the system reaches a clustered state. Otherwise, the system is considered as in a gaseous phase.

Figure 2 gives examples of how the K–S test method is applied to find the cluster state in our simulation observations. Three examples are given in Figs. 2(a)–2(c) with the cell length and the total number of particles: Lx = 75 and N = 3000, Lx = 80 and N = 5000, and Lx = 30 and N = 4000, respectively. In each case, the particle density distribution function ρ(x) (left), its cumulative distribution function F(x) (middle), and the front view of cell with particles’ Voronoi volume fraction (ϕlocal = Vp/Vvoro) (right) are plotted in Fig. 2. In this paper, the Voronoi volume fraction Vvoro is calculated by Voro++ code.[18] When the count of particles with ϕlocal ⩾ 0.285 appears and has a sudden increase, the system is considered as in cluster phase. These two criteria agree with each other when the cluster width is smaller than δ/2. For dense system, when the cell length δ/2 is smaller than the cluster size, K–S test does not apply and we consider the system in a dense cluster state. For example, in the cases of Figs. 2(a)–2(c), Da = 0.130, Db = 0.192, and Dc = 0.154. For α = 0.05, Kα = 1.358. In Fig. 2(a), we have n = 75, , , and particles’ ϕlocal are smaller than 0.285. It is thus a gas state. For the case in Fig. 2(b), n = 80, , , and ϕlocal ≥ 0.285 particles appear, thus it is considered as a cluster state. For Fig. 2(c), n = 30, , , but with a large number of ϕlocal ⩾ 0.285 particles. It is taken as in a dense cluster state.

Fig. 2. Left: the particle density distribution ρ(x) (black scatter). Middle: the cumulative distribution F(x) (black solid line), and the cumulative distribution function of a uniform function U(x) (red dotted line). Right: the front view of the cell with particles’ Voronoi volume fraction. For different simulation parameters Lx and N, (a) gas phase, (b) cluster phases, and (c) dense cluster are shown. (a) Lx = 75, N = 3000, (b) Lx = 80, N = 5000, (c) Lx = 30, N = 4000.
3. Results
3.1. Heating in gas state

All the analyses given above are based on the system reaching a steady state in the granular gas phase. We then need to know the time for the system to reach a steady state. Here, in Fig. 3(a), we show a typical example of global temperature versus time using parameters in the reference set in Table 2. A stretched exponential function[19]

is used to fit with the growth T(t) (see Fig. 3(a)). It shows that a typical time τα for the system to reach a gas steady state is 1.8 s in the case of Fig. 3(a). To estimate the influence of δ/r and ϕg on this heating time, the characteristic heating time τα is investigated and the result is given in Fig. 3(b). For all the cases in the reference set (Table 2), fitting parameter c = 2, and the time τα changes in the range of 1.8–2.8 s. As δ/r increases, the energy transport slows down as τα increases to ∼ 2.8 s. The global volume fraction ϕg has only a slight effect on τα. Increasing the volume fraction ϕg would decrease the heating time τα. In Fig. 3(b), the largest τα (4 s) appears at ϕg = 0.0068. To make sure the system reaches the steady state, all simulation data points are taken at a time t ⩾ τα, normally after 6 s.

Fig. 3. (a) Example of global temperature T (set No. 1 in Table 2) versus time, and its fitting by a stretched exponential function (Eq. (8)) (red line). (b) A map of the fitting parameter τα for those gas state data points in the δ/r–ϕg phase diagram with the reference set (set No. 1 in Table 2).

In this work, for each set of parameters in Table 2, we run at least 55 data points. Of the seven sets in Table 2, the first set of parameters is the same as those in microgravity experiments and is used as a reference set.

3.2. Cooling to cluster state

In this simulation, we only consider cluster forming due to inelastic collisions among particles. In this section, we investigate the free cooling[20,21] process without energy input and with no sidewall collisions (periodic boundary condition as set No. 7 in Table 2). The global granular temperature is calculated as the averaged fluctuations of particle velocities in three dimensions

The initial temperature of the system is set to be T0 = T(t = 0). The temperature T(t) decreases rapidly over time through particle collisions as shown in Fig. 4(a). The decay is well fitted by the Haff’s law[22] , where the fitting time τH = 6.58 s in Fig. 4(a) when the restitution coefficient ε = 0.1. In Fig. 4(b), the characteristic cooling time is plotted as a function of ε. We have modified Haff’s time by adding a friction term γ to make the time τH given by[6]

where v0, ε, η = N/V, σ = π (2r)2, and C0 are the initial velocity, the restitution coefficient of the particles, the number density of the system, the cross-section of a particle, and a calibration parameter, respectively. The fitting result shown in Fig. 4(b) gives γ = 0.24.

Fig. 4. (a) Simulation results of temperature decay with time and the Haff’s law fitting (red solid line) with the fitting parameter τH = 5.68 s. (b) Typical cooling time (τH) versus restitution coefficient, with fitting parameters C0 = 4.82, γ = 0.24 (red solid line).
3.3. Phase diagram

Considering friction coefficient 0.1 and restitution coefficient 0.9, simulation results are plotted in a δ/r–ϕg diagram, as shown in Fig. 5. A simple model[16] based on the competition between the cooling time τH and a propagation time τp is used to the transition curve of the gas–cluster phase diagram. The time τH is given in Eq. (9).

Fig. 5. The phase diagram δ/r versus ϕg using the first set of parameters in Table 2. Triangles are in the gas state, and circles are in the cluster state. The red solid fitting curve is a best fit with Eq. (11), in which the restitution coefficient ε = 0.9 and the fitting parameter C0 = 10.5 (red line). The black and blue dashed lines give a range of fitting for C0 = 11 (black dashed line), C0 = 10 (blue dashed line). The three points a, b, and c are the feature points corresponds to Figs. 2(a)–2(c), respectively.

The propagation time τp can be calculated by simplifying the problem as a one-dimensional granular chain with equal intervals Δx

which is a time for a particle with an initial velocity v0 to transmit from one cell end to the other end, that is, the distance δ between the two pistons.

The propagation time τp (Eq. (10)) is equal to .[14] When τH ≤ τp, particle may lose its velocity before reaching the other piston, and cluster forms. The competition between the cooling time and the propagation time provides us the condition of gas–cluster transition, that is τH = τp,

where r is the radius of the particle 0.5 mm, ϕg is the global volume fraction, R0 = 2, and γ = 0.24. Here, C0 is a calibration parameter.

Fitting the data by Eq. (11) will yield the gas-cluster transition curve with C0 as a fitting parameter. We find C0 = 10.5 is best fitted as shown in Fig. 5. The dense cluster region is colored in purple in the δ/r–ϕg diagram of Fig. 5. In the following, we will investigate the parametric influence on the transition phase diagram.

3.4. Parametric influences

Using K–S test and Voronoi volume fraction, we have analyzed the simulation data and obtained the gas-cluster phase diagram in δ/r–ϕg plot. Simulations with periodic and fixed sidewall boundary conditions are performed. The friction coefficient μpp = 0.1, μpw = 0.1, and the restitution coefficient ε = 0.9 are chosen. In Fig. 6, comparison of cluster distribution under PBC and fixed boundary condition for a cell with Lx = 70 mm, N = 5000 is given. The slice shown in Fig. 6 is cut at x in between 30 mm and 40 mm. The chromaticity diagram shows ϕlocal in the range of 0.05–0.65. It is shown that for PBC, cluster forms in the central area, while in fixed sidewall boundary condition, the clusters tend to initiate from the corner or along the sidewalls, as seen in Fig. 6(b). In order to eliminate the boundary effect, we choose PBC for the remaining parametric studies.

Fig. 6. The comparison of cluster distribution for PBC and fixed boundary condition. The simulation parameters are μpp = μpw = 0.1, ε = 0.9, and Lx = 70 mm, N = 5000. The slice shown in the figure is chosen in the middle of the cell, Lx = 30 mm–40 mm. The chromaticity diagram shows ϕlocal from 0.05 to 0.65.

The influence of three material parameters, particle–particle friction coefficient μpp, particle–wall friction coefficient μpw, and the restitution coefficient ε, are investigated in the following simulations. In Fig. 7(a), as the particle–particle friction coefficient μpp increases from 0.1 to 0.5, the fitting parameter C0 changes to 8.5 from 10.5, i.e., the gas–cluster transition shifts to smaller ϕg. When changing particle–wall friction coefficient μpw from 0.1 to 0.5, the influence is nearly negligible (as shown in Fig. 7(b)). Figures 7(c) and 7(d) show the influence of the restitution coefficient ε. It is seen that for ε = 0.97, the value of fitting parameter changes to C0 = 5 from C0 = 10.5 when ε = 0.9. The transition curve shifts to greater ϕg, that is, a denser system with more collisions is needed for the cluster to occur. For ε = 0.1, the enormous C0 indicates that the fitting fails in this extreme case. In addition, a phase diagram by a fixed boundary simulation is shown in Fig. 7(e) to compare with the simulation results by periodic boundary phase diagram shown in Fig. 5. The fitting curve shifts from C0 = 10.5 to C0 = 9. This indicates that there exists some minor sidewall boundary effect.

Fig. 7. (a) Phase diagram with parameters in set No. 2 of Table 2. (b) Phase diagram with parameters in set No. 3 of Table 2. (c) Phase diagram with parameters in set No. 4 of Table 2. (d) Phase diagram with parameters in set No. 5 of Table 2. (e) Phase diagram with parameters in set No. 6 of Table 2. The green area is gas state, and the red area is cluster state. The triangles are in the gas state, and the circles are in the cluster state. The red solid line is the fitting curve using Eq. (11), γ = 0.24, R0 = 2, and C0 is fitting parameter.
6. Conclusion

In this work, we have performed a simulation study on the parametric influences of the clustering transition in a vibration-driven granular gas system. The parameters studied include the global volume fraction of the system, the size of the system, the friction coefficient, and the restitution coefficient among particles and among particle–walls. Periodic boundary and fixed boundary of sidewalls are also checked in the simulation. It is found that in our model, the restitution coefficient of the particle-particle collision plays the major role in influencing the transition. Friction among particles plays a minor role, and a friction term γ is needed to modify the Haff’s cooling time. This simulation study helps us in the determination of the experimental parameters in the future microgravity experimental preparations.

Reference
[1] Jaeger H M Nagel S R Behringer R P 1996 Rev. Mod. Phys. 68 1259
[2] Kadanoff L P 1999 Rev. Mod. Phys. 71 435
[3] Gennes P G D 1999 Rev. Mod. Phys. 71 S374
[4] Hou M Liu R Zhai G Sun Z Lu K Garrabos Y Evesque P 2008 Microgravity Sci. Technol. 20 73
[5] Hou M Li Y C Liu R Zhang Y Lu K Q 2010 Phys. Status Solidi (a) 207 2739
[6] Wang W G 2018 Chin. Phys. 27 084501
[7] Miller S Luding S 2004 Phys. Rev. 69 031305
[8] Falcon E Wunenburger R Évesque P Fauve S Chabot C Garrabos Y Beysens D 1999 Phys. Rev. Lett. 83 440
[9] Pathak S N Jabeen Z Das D Rajesh R 2014 Phys. Rev. Lett. 112 038001
[10] Falcon E Fauve S Laroche C 1999 Eur. Phys. J. 9 183
[11] Kudrolli A Wolpert M Gollub J P 1997 Phys. Rev. Lett. 78 1383
[12] Noirhomme M Ludewig F Vandewalle N Opsome E 2017 Phys. Rev. 95 022905
[13] Opsomer E Ludewig F Vandewalle N 2011 Phys. Rev. 84 051306
[14] Opsomer E Ludewig F Vandewalle N 2012 Europhys. Lett. 99 40001
[15] Kloss C Goniva C 2010 Proc. 5th International Conference on Discrete Element Methods August 25--26, 2010 London, UK
[16] Noirhomme M Cazaubiel A Darras A Falcon E Fischer D Garrabos Y Lecoutre-Chabot C Merminod S Opsomer E Palencia F Schockmel J Stannarius R Vandewalle N 2018 Europhys. Lett. 123 14003
[17] Sachs L 1997 Einführung in die Statistik. In: Angewandte Statistik Berlin Springer 427 10.1007/978-3-662-05746-9_2
[18] Rycroft C 2009 Voro ++: A three-dimensional Voronoi cell library in C++ Lawrence Berkeley National Lab. (LBNL) Berkeley, CA 10.1002/pssa.201026461
[19] Ciamarra M P Coniglio A Nicodemi M 2006 Phys. Rev. Lett. 97 158001
[20] Maaß C C Isert N Maret G Aegerte C M 2008 Phys. Rev. Lett. 100 248001
[21] Harth K Trittel T Wegner S Stannarius R 2018 Phys. Rev. Lett. 120 214301
[22] Haff P K 1983 J. Fluid Mech. 134 401